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We study the simulation of charged systems in the presence of general boundary conditions in a 
local Monte Carlo algorithm based on a constrained electric field. We firstly show how to implement 
constant-potential, Dirichlet, boundary conditions by introducing extra Monte Carlo moves to the 
algorithm. Secondly, we show the interest of the algorithm for studying systems which require 
anisotropic electrostatic boundary conditions for simulating planar geometries such as membranes. 



I. INTRODUCTION 

In the simulation of condensed matter one very often 
imposes periodic boundary conditions in order to mini- 
mize surface artefacts which give rise to slow convergence 
of energies to their bulk values. While such boundary 
conditions are simple to understand for short-ranged po- 
tentials they lead to many subtleties when working with 
charged systems^— . Attempts at simplifying the prob- 
lem by using truncated potentials lead to violations of 
basic sum rules on the structure factor—, incorrect num- 
ber fluctuations in finite systems^ or even loss of electrical 
conductivity^. The correct treatment of charged media 
requires a full treatment of the long-ranged Coulomb in- 
teraction. 

Classical methods of treating the Coulomb interaction 
use careful mathematical analysis to convert the con- 
ditionally convergent Coulomb sum into a well defined 
mathematical object such as the Ewald sum 1 -. Recently 
we introduced an alternative treatment of the Coulomb 
problem that allows one to replace the global calculation 
of the interaction energy by a local dynamic process&i^. 
The main interest in this transformation is that it enables 
one to perform a local Monte Carlo simulation in the 
presence of charges and dielectrics, without ever solving 
the Poisson equation. The algorithm works by evolving 
the electric field in time (in a manner similar to Maxwell's 
equations) while eliminating "uninteresting" degrees of 
freedom such as the magnetic field. The algorithm uses 
the energy 
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(1) 



where E is the electric field while implementing Gauss' 
law 



e V • E - p = 



(2) 



as a dynamic constraint. We showed^ that the algorithm 
can be used to generate tin-foil or vacuum boundary con- 
ditions if one chooses appropriate dynamics for the q = 
component of the electric field. Until now applications 



have been to the properties of bulk, three dimensional 
medi a 10 ' 11 in periodic systems. 

This paper generalizes the method to a broader class 
of physical systems and boundary conditions. Firstly we 
consider the case, particularly important in devices and 
electrodes, of imposition of an external potential on a 
metallic surface; a case which requires the use of Dirichlet 
boundary conditions for the potential. Again we find that 
the locality of the formulation allows the simulation of a 
broader class of geometries, including those for which the 
use of fast Fourier techniques is difficult. 

A second generalization of the method is required 
to treat anisotropic systems, in particular membranes. 
The simulation of thin, quasi two-dimensional systems 
in three-dimensional space is surprisingly difficul t 12 ' 13 ' 14 . 
The use of a purely local algorithm requires only small 
modifications in order to minimize finite size effects; we 
argue that it will also give rather favorable complexity in 
the limit of large number of particles, N. 

Our paper contains two, independent, self-contained 
theoretical sections. Firstly in II we treat the problem 
of metallic boundary conditions including certain sub- 
tleties as to how a true metal behaves. Secondly in III 
we consider electrostatics in anisotropic 2 + 1 dimensional 
geometries. Implementation of the ideas has been per- 
formed by Thompson and Rottler— using off lattice tech- 
niques suitable for atomistic simulation. In their paper 
they compare detailed simulations with available theo- 



II. METALLIC BOUNDARY CONDITIONS 

We begin by emphasizing the difference between ide- 
alized metallic boundaries and true physical systems in 
which screening occurs over a small but finite Debye 
length. We then show how to impose Dirichlet condi- 
tions on the potential, 4> while still keeping E as the main 
dynamic variable. 
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The nature of metallic boundary conditions 

There are two different, seif-consistent ways of calcu- 
lating the energy of charges in the presence of a dielectric 
or conducting interfaces. In the first treatment we calcu- 
late the potential energy arising from the solution of the 
Poisson equation 

V • (eV0 = -p 

with p the density of free charges, with the appropriate 
boundary conditions for each configuration. The bound- 
ary condition at a metallic surface is that the tangential 
electric field vanishes, while the normal field at the sur- 
face satisfies 



D n = eE n = — a 



(3) 



where a is the surface charge density. The field is iden- 
tically zero within a conductor—. The internal energy is 
then 



U 



>d 3 7 



(4) 



In this case the longitudinal electric field is static, and 
is given by E = — V(/>. We note that in this treatment 
thermal fluctuations play no role so that certain fluctua- 
tion phenomena such as thermal Casimir interactions are 
neglected. 

A physical system at a non-zero temperature always 
presents polarization or charge fluctuations; superposed 
on the interaction energy Eq. Q is an effective poten- 
tial coming from these fluctuations. We now work in the 
Debye-Huckel limit to calculate an approximate correla- 
tion function for the electric field in a conducting sys- 
tem at finite temperatures. This illustrative calculation 
shows that transverse field fluctuations become impor- 
tant as soon as the temperature T > 0. 

We consider a one component plasma, with neutral- 
izing background. We approximate the free energy of a 
conductor as a sum of the electric field energy Eq. (TT]) 
and the configurational entropy of the charges with the 
functional^ 



F = 



UbT cine > d r 



(5) 



where c(r) = p(r)/e is the number density of mobile ions 
and e the charge of the particles. In the limit of small 
fluctuations of density Sc(r) we expand to second order 
and use charge conservation so that / (5cd 3 r = 0. Then 



F(<5c,E) 



2c 
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where Co is the mean background density of the ions. 
We now eliminate the charge fluctuation with the help 
of Gauss' law and find an effective action for the electric 
field 



Fe = e i — 



E 2 
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with k the inverse Dcbye length. From Eq. ([6]) we read off 
the longitudinal and transverse correlations of the electric 
field, 

k T 

e„{E{q)-E(q)) long = y^y^, 
eo(E(«)E(g)) tran — kgT. 

On wavelengths large compared with the screening length 
(q/ K « 1) all longitudinal and transverse electric modes 
fluctuate with the same amplitude. 

It is instructive to compare with similar calculations 
for an explicit model of dielectric in terms of polarization 
fluctuations^. We find that 



e (E(g)E(g)) 
e (E(q)E(q)) t 



long 



k B T 
= k R T. 



Again transverse field correlations are unchanged in the 
presence of a dielectric or classical charged fluid, longitu- 
dinal fluctuations depend on the material properties. In 
the limit e — > oo the fluctuations of a dielectric are the 
same as those of a metal for q <C K. 

From these expressions in Fourier space we calculate 
the correlations of the electric field in real space. For a 
dielectric (with e < oo) we find that the field displays 
dipolar correlations so that 
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field correlations are long-ranged. This dipolar corre- 
lation in the fields leads to long-ranged Casimir /Lifshitz 
interactions between dielectrics. For a conducting system 
field correlations decay exponentially with characteristic 
length the Debye length, 

( J B i (r) J B j (r')}~e-' t l r - r 'l. 

From these considerations we sec that different strate- 
gies are needed to simulate the idealized classical metallic 
boundary condition Eq. ([3]) or a true fluctuating charged 
fluid. In the first situation one only needs information 
on the fields in the nonconducting regions. In the sec- 
ond case the surface field is influenced by fluctuations 
that occur within a few Debye lengths of the surface. 
The surface must then be simulated explicitly with an 
expression such as Eq. §5§ in order to obtain results in- 
cluding Casimir/Lifshitz type interactions. We note that 
such interactions are very weak on the macroscopic scale. 
Between two plates of area A separated by a distance H 
the thermal interactions are given bj*i£ 
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C(3) 
16ttH 2 



with C(3) « 1.20 

In this paper we will only consider the case of imposing 
the first type of metallic boundary conditions in which 
fluctuation forces are neglected. 
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Fixed-potential boundary conditions 

Minimization principle 

We now apply our local simulation algorithm to fixed- 
potential boundary conditions of the type Eq. © . The 
method of attack is a generalization of the methods of 
Ref. 0. Firstly we seek a variational principle for the elec- 
tric field for which the energy is a true minimum for the 
geometry of interest. We then promote the constraints of 
the minimization principle to (^-functions in a partition 
function. This partition function is then shown to gen- 
erate the same relative statistical weights as the original 
variational energy, but is much more convenient for the 
purposes of simulation since it allows the local update of 
fields in a Monte Carlo algorithm. 

As noted in the introduction electrostatic interaction is 
calculated by minimizing the electric field energy Eq. ((T|) 
in the presence of the constraint of Gauss' law, Eq. (2j, 
using the functional 



A[E] 



e E 2 



A(e V ■ E - p) 



d 6 r, 



where A is a Lagrange multiplier. We now generalize 
Ref. Qj^ and consider a system with both free charges 
and conductors, i. These conductors are maintained by 
external sources at fixed potentials </> l -° xt ' ) . This is most 
conveniently done by performing Legendre transform^ in 
order to eliminate for the unknown surface charge densi- 
ties a = {<Ji} on the surfaces S = {Si} in order to replace 
them by the potentials 0( cxt ) = {(f>[°^}- 



C/ 2 [E,(t] = U[E] - & a^^dS. 



(7) 



We now impose the boundary condition Eq. ([3]) with a 
further Lagrange multiplier p living on each surface and 
consider the functional 



A 2 = U 2 



J A(e V-E-p)d'V+^V(e n-E + cr) dS. (8) 

In order to find the stationary point of A 2 we vary all the 
fields, and integrate by parts using AV ■ <5E = V • (A(5E) — 
(VA) ■ 8E and find 

SA 2 = J [e £E • (E — VA) + SX (e V ■ E - p)] d 3 r 
- <£[5<j{<j) {cxt) +p)+5p(e n-E+<T)-e {\-p)SE-n] dS. 
At the stationary point we find 
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The first equation implies that E is the gradient of a 
potential, 4> = —A. The next two imply that <p = 0( cxt ) 
on S so that we indeed generate Dirichlet conditions. 
Finally the last two equations impose Gauss' law. At the 
stationary point thus we have E = E p = — V^> p and a = 
a p , where the index p signals the solution to Poisson's 
equation with Dirichlet boundary conditions. 



Partition function 

We now generalize the stationary principle to finite 
temperatures and replace the Lagrange multipliers of 
Eq. (|8|) by ^-functions in a functional integral. Let us 
define a partial partition function, where the bulk charge 
distribution p is given, as 

Z p = J VEV<j <5(e V-E-p) ( 5(eon-E+CT)e-' 3C/2[E ' <Tl . (9) 

We integrate over the electric field, constrained by Gauss' 
law and also all values of the surface charge a compatible 
with the flux condition at the surface of the conductors. 

The full partition function is calculated by integrating 
over the position of all mobile charges. In order to find 
the interaction in terms of the minimum of the functional 
it is convenient to change integration variables so that 
E = E p + e and a = a p + s. After noting that a p = 
—eon • E p we find 



Z p = Jve Vs <5(V • e) <J(e n ■ e + s) 

x e -P[S l(E P +e) 2 d 3 r- f(a p + s)^ xt) dS] 



With the help of the constraint equations we find 

-PU 2 [E P ,<J P ] 

,-PU 2 [e,s] 



Z P = c 



J Ve Vs <5(V • e) S(e n • e + s) 



so that 



Z p = C -PU 2 IE P ,* P ] z 



fiuct i 



(10) 



with Zfl uc t independent of positions of the mobile charges. 
This results in a complete factorization of the statistical 
weights of charge configurations and of field and surface 
charge fluctuations 



Z = 



-/3(7 2 [E p ,cr p ] 



^fluct 



^Coulomb-^fiuct 



We note, however, that following the discussion of Section 
II, this fluctuation partition function is not necessarily 
the true physical fluctuation partition function that one 
would calculate for a true physical metal. Thus sampling 
with the partition function Eq. ([9]) will not necessarily 
give the correct fluctuation forces acting on electrodes 
even if, by construction, it generates the correct relative 
weight for configurations of free charges. 
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FIG. 1: Local update for integrating surface charge fluctua- 
tions. Grey region represents the conductor volume (where 
E = 0); the white region is the dielectric. 



Algorithm 




FIG. 2: Updating the total charge of conductors implies up- 
dating the field along a line connecting them. This update is 
not local. 



The result Eq. (|10|) is the analytical basis for an algo- 
rithm that simulates fixed-potential surfaces. Sampling 
the partition function requires updates generating fluctu- 
ations of the variables, {ij}, E and a. We refer the reader 
to previous work&&2ii for particle updates and bulk field 
updates. We note that the demonstration requires an 
integral over the surface charge er; thus even if charges 
in the volume are discrete those on surfaces should be 
sampled continuously. 

In order to sample the partition function Eq. ^ two 
new moves need to be introduced. Firstly, fluctuation in 
the charge density in the surface of a given conductor that 
conserve its total charge. This move is implemented by 
a variation of the simplest update for volume charges. A 
random charge amplitude 5 is generated to make a charge 
pair of +5 at one site and —S at a second site. Since field 
lines must lie outside the conductor, Fig. [H we update 
the field on three links connecting the modified sites. To 
sample the integral Eq. @ the charges should be chosen 
from a continuum distribution. Secondly, each conductor 
i bears a net charge J s <ji which should also fluctuate. 
This is achieved by transferring charges between pairs of 
conductors. Gauss' law then requires modifying the total 
electric flux between the two surfaces. 

We now demonstrate that despite the need to trans- 
fer charge large distances between conductors this does 
not dominate the time need to simulate the system. We 
start by noting that fluctuations on a conductor can be 
estimated using an argument based on the capacitance 
of an isolated conductor—. In three dimensions the ca- 
pacitance of an object of size d scales as C = e d. If we 
equate the charging energy to the thermal energy scale 
we find 
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e d 



so that the amplitude of charge fluctuation of a conduc- 
tor is Q 2 ~ ksTeod. Consider a Monte Carlo trial in 
which a charge 5 is sent along a path of length £, Fig. [2] 
Ther>2£ we estimate the energy change as £S 2 /eoa 2 , with 
a the mesh spacing. For this trial to be successful we 



again require that this energy is matched to the thermal 
fluctuations giving 8 2 ~ eoa 2 kBT/£. Thus the amplitude 
of the charge transfer over large distances must be small. 

We now consider M such exchanges. Since charges are 
transferred in each direction we require that MS 2 = Q 2 , 
or M ~ Id I a 2 . Each transfer requires a numerical effort 
that is also 0(£/a) to update the links so that a total 
computational effort T) = £ 2 d/a 3 is needed to equilibrate 
the charge fluctuations between the conductors. In any 
practical simulation in a box of dimension L we expect 
that n < L 3 /a 3 , so that the effort need to equilibrate 
charges between conductors is less than that required to 
perform a single sweep of the bulk of the simulation. 

We thus have the basis for the generalization of a lo- 
cal Monte Carlo algorithm for simulating charges in the 
presence of Dirichlet boundary conditions. Conventional 
methods for treating this problem use fast Fourier tech- 
niques in simple separable geometries, such as regular 
simulation cells. The use of a local formulation of the 
electrostatics allows one to implement such boundary 
conditions in arbitrary geometries — including rough sur- 
faces or irregular electrodes where Fourier analysis does 
not work. 



III. SIMULATING SYSTEMS WITH SLAB 
GEOMETRIES 

We now pass to consideration of a second important 
class of boundary conditions that should be imposed in 
the study of systems with planar geometries. Exam- 
ples include thin quasi two-dimensional slabs isolated in 
three-dimensional space, or in the biophysical field sim- 
ulation of proteins in a lipid membrane with an implicit 
solvent. In order to generate the correct Coulomb in- 
teractions between particles with grid based methods 
the thin sample must be embedded in a thick three- 
dimensional simulation box. This box then contains 
many more degrees of freedom than the original parti- 
cle system. For reasons of efficiency one wishes to make 
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the repeat distance perpendicular to the thin sample as 
small as possible, however if this repeat distance is too 
small multiple copies of the sheet "see" each other and 
introduce artefacts in the simulation. 

Imposition of tin-foil or vacuum boundary conditions 
gives rise to simulations that converge very slowly as the 
system size is increased. The crucial insight into how to 
accelerate this convergence was provided in Ref. [lU and 
efficient implementations are now available that translate 
these idea o 14 i 22 in molecular dynamics. Our aim here is 
to reproduce this rapid convergence in a local formulation 
of electrostatics, and to argue that the asymptotic com- 
plexity is at least as good as TV 3 / 2 and in many systems 
N 1 . 

In standard approaches to the electrostatic energy all 
the complexity comes from the transformation from an 
ambiguous, conditionally convergent sum (a sum whose 
answer depends on the order of evaluation) into a rapidly 
converging, unambiguous expression, such as the Ewald 
formulai. In local formulations there is no equivalent to 
the conditional convergence. Instead one is left with dif- 
ferent, inequivalent, choices for the treatment of the zero 
wavevector component of the electric field ~E(q = 0)2. Let 
us first review the origin of these different choices. 



Periodic boundary conditions in local algorithms 

In periodic boundary conditions an arbitrary vector 
field can be decomposed into three terms 

E = + V x G + E 

where the potential (f> is periodic, as is G. E is a constant; 
it is the zero wavevector component of the electric field, 
E(g = 0). With this decomposition the energy is given 
by the sum of three independent terms 



U 
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V E^ 



with V the volume of the simulation cell. Cross-terms 
are shown to be zero on integration by parts. The sim- 
plest version of the local Monte Carlo algorithm uses two 
updates^. 

• Motion of charges q along links I together with a 
slaved update of the electric field on the link ac- 
cording to Ei — ► Ei — q/eo- Both the longitudinal 
and transverse components of the field are modi- 
fied, as well as the q = component, E. 

• A plaquette update which leaves unchanged both 
E and the longitudinal field — V0, modifying only 
G. 

This pair of updates are respectively a discretized analog 
of the two terms on the right hand side of the Maxwell 
equation 



eo 



<9E 

dt 



Integrating the Maxwell equation over both time and the 
simulation volume we see that the q = component of the 
field and the dipole moment of the system of charges d = 
Jdt Jd 3 rJ are linked by eo^E + d = const, a conservation 
law which remains valid for the discretized equations. 

Thus the simplest form of the algorithm samples the 
partition function 

Z= [VE < 5(V-E- /9 )(5(eoFE + d)e-' 3 -'' eoE2/2d3r . (11) 



If we introduce the solution to the Poisson equation in 
periodic boundary conditions <p p then we see that this 
partition function samples configurations with the effec- 
tive energy 



Ueff = 



2e V 



(12) 



This result is very similar to that found using careful 
evaluation of the Coulomb sumi which gives 



Ulps 



2(l + 2e s )e V 



where p is the cell dipole moment and e s the relative 
dielectric constant of an "exterior" reference medium, 
two common choices being e s = 1 for vacuum boundary 
conditions and e s = oo for metallic or tin-foil boundary 
conditions. The expression resulting from the local al- 
gorithm is very similar to that found in the summation 
method if we take e s = and we notice that d and p are 
almost identical, p is the dipole moment of the charge 
images inside the Bravais cell being used in the simula- 
tion; when a particle crosses the boundary of the cell it 
is replaced by that of its periodic images which enters 
the cell so that p is discontinuous, d is a dipole moment 
which is continuous on crossing Bravais cell boundaries, 
thus it keeps track of how many times the particles wind 
around the simulation box: 



J + V x H. 



where qi is the charge of each particle in the simulation 
cell and is a Bravais lattice vector. 

The authors of Ref. HI introduced another update 
scheme for the electric field based on a worm algorithm 
widely used in the simulation of quantum spins^ 3 -. The 
algorithm nucleates a pair of virtual charges ±q v . One 
of these charges diffuses until it returns to its companion 
where the two charges annihilate. If the mobile charge is 
confined to the periodic cell then during the move only 
the transverse field V x G is updated. If the particle is 
allowed to wind around the cell it breaks the conserva- 
tion law on eoT^E + d, removing the second (5-function in 
Eq. (fTT|). We then find that the effective energy Eq. (TJ) 
is independent of d; we have effectively periodic, that is 
tin-foil, boundary conditions. A similar result is found 
by including E as a single extra dynamic variable which 
must be updated with a third, independent Monte Carlo 
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Finite size effects in slab geometries 

We now give a short argument for the slow convergence 
of the interaction energy for systems which are simu- 
lated in slab geometries. Rigorous calculations which 
provide the full justification are to be found in the 
literatur o 2 ' 12 i 14 . Consider a thin, quasi two-dimensional 
system embedded in a simulation box of dimensions 
L 2 x L z ; usually one is interested in the case L z /L > 1 
so that successive periodic images do not interact too 
strongly. 

The slow convergence of the energy can be understood 
by considering the electrostatic potential in a mixed real- 
space/Fourier representation, (f>(q\\,z) where denotes 
the transverse wavevector in the (x, y) plane. Using the 
fact that polarization of a sample is equivalent to a space 
charge p = — V ■ P the mean polarization of the thin sam- 
ple in the z direction is equivalent to two charge sheets 
of strength a = ±p z /L 2 h where p z is the z component 
of the dipole moment and h the sample thickness. The 
interesting physics occurs in the equation for = 0. We 
consider a single cell in the z direction, placing the source 
—a at z = and the second +a at z = h. Then 

d 2 cf>(0,z) 
£ ° — ~dz~ 2 — = ° ~ ^ ~ ^ ' 

This has as a solution 

0(0, z) = (-e + cr/e )z 0<z<h 
= {— eo + c/e )/i - e (z — h) h < z < L z 

where we have used continuity, but not periodicity, of 
the potential, eo is a yet to be determined integration 
constant; it corresponds to the electric field flowing be- 
tween two copies of the simulation cell, Fig. [31 The z 
component of the electric field is 

E z = E mt = e - cr/e < z < h 

= E ext = e h< z <L. 

The energy of this simulation cell is the integral of 
e E 2 /2, 




The classical tin-foil solution is recovered by setting eo = 
ah/eoL z so that the potential is periodic, </>(0) = <p(L z ), 
Fig. [3j and the average field E z = J dz E z (z) = 0. 
But the energy then depends on L z , 

U Coulomb = Uoo — - — Pz (13) 

2e L z L z 

where Uoo is the energy in the limit of large L z . If one 
considers a series of slabs with identical dipole moment 
and transverse dimensions and vary L z then the energy 
converges as only 1/L Z ; very large, and mostly empty 
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FIG. 3: Plot of (p{0, z) in the presence of a double layer. 
E mt = eo — cr/eo, E ext = eo, with eo = ah/eoL z . Poten- 
tial plotted over several simulation cells. The non-zero E ex 
couples successive periodic images. 
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FIG. 4: Plot of $ with eo = 0. The field between slabs, E ext , 
is zero. There is no dipolar interaction between the slabs. 



boxes must be simulated in order to reduce finite size 
artefacts. 

There is no such convergence problem when using 
anisotropic electrostatic boundary conditions in which 
the exterior field is always zero: eo = 0, Fig. 2J The 
energy is Ucoulomb — U^, independent of L z . The local 
algorithm automatically gives eo = 0, as do Maxwell's 
equations; in order to still have E x = E y = we re- 
quire tin-foil boundary conditions in the (x, y) plane. Our 
method implicitly includes the +p 2 /2tQ,V correction that 
potential-based algorithms need to add to the tin-foil en- 
ergy Since the particles remain confined to a single lat- 
tice cell in the z direction there is no distinction between 
d z and p z . 

Let us now consider the effect of density fluctuations 
with the plane, by studying the Poisson equation for 
wavevectors qp ^ 0: 

2 ~ d 2 $ 
0+ c^ = 1— 
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FIG. 5: Electrostatic energy of the slab system when all three 
components of E are sampled independently, corresponding 
to tin-foil boundary conditions. Simulation results (+) follow 
Eq. ([15) . L = 15, e = 1. Five values of k B T from 0.1 to 0.5 
for each L z . 

Outside of the slab where p — we conclude that 

0(q,l ^0,z> h) ~e ±q " z 

where q|| = corresponds to the special case treated 
above. The next longest-ranged component to the inter- 
action must come from modes of the form (27r/L,0,0) 
and lead to interactions decaying as e _27rLz / L . Already 
when L z /L = 3 one finds a reduction in the interaction 
by a factor 7 x 10~ 9 . 
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FIG. 6: Coulomb energy of the slab system when E z is not 
updated independently of charges. Simulation results (+) arc 
independent of L z . L = 15, eo = 1. Five values of ksT from 
0.1 to 0.5 for each L z . 

energy varying as 1 /L z (with the analytically determined 
prefactor). In a second series of simulations updates to 
E z are dropped. Fig. [6] shows that the electrostatic en- 
ergy is now independent of L z to within the precision of 
the measurements and equal to Uoq. The natural cou- 
pling between d and E introduced by our local method 
removes the dipolar interaction between replicas of the 
slab. 



Discussion 



Numerical tests 

To test these ideas we simulated a lattice gas with 
a slab geometry using the local algorithm. The lattice 
spacing a is set to unity. We start a simulation with all 
positive and negative charges superposed at z = and 
E = in the whole simulation volume. We then displaced 
positive charges to z = 1, creating a dipolar sheet, up- 
dating the fields according to the constrained algorithm. 
We simulated the fields using several different tempera- 
tures, T and cell heights, L z . The total energy can be 
decomposed into a static, electrostatic contribution plus 
a thermal energy N p iksT '/2 coming from equipartition in 
the N p i degrees of freedom associated with the transverse 
excitations. On fitting the total energy to the form 

(U) = Ucoulornb{L z ) + — - P ' 

we extracted an estimate of the electrostatic energy as a 
function of the cell dimensions. 

In a first series of simulations we imposed the equiva- 
lent of tin-foil boundary conditions by including the three 
components of E in the set of dynamic fields. The results 
are shown in Fig. [51 as a guide to the eye we added the 
continuous curve which corresponds to a correction in the 



In a local Monte Carlo simulation the most obvious 
implementation would use a mixture of local updates of 
the particles together with the worm algorithm-'- 3 in or- 
der to integrate over the transverse degrees of freedom in 
the empty bulk of the simulation. This worm algorithm 
re-equilibrates the transverse field with an effort which is 
O(l) per updated degree of freedom (which corresponds 
to plaquettes). If one were to use just a single worm up- 
date per sweep of the particles then we find a complexity 
per sweep which varies as 0{N) for the particle motion, 
and 0(N 3 / 2 ) for the worm dynamics. In practice launch- 
ing a worm which simulates the entire simulation box for 
every particle sweep results in an oversampling of the 
uninteresting transverse degrees of freedom. If one simu- 
lates a set of particles for a time r then diffusion motion 
in the plane will give a typical displacement l t ~ \Jt. A 
worm which updates plaquettes a distance more than t t 
from the slab is doing unnecessary work. 

This suggests the following hierarchical scheme: every 
sweep we launch a worm confined to the simulation slab 
thickness h, then every 4 sweeps we launch a worm con- 
fined to a distance 2h. Similarly every 2 2n sweeps we 
allow the worm to propagate 2 n h in the z direction. In 
such a scheme the longest wavelength modes of the elec- 
tric field are re-equilibrated on the same timescale as is 
needed for particles to diffuse across the simulation box. 
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The worm then spends most of its time updating at the 
scale of the slab so that the total algorithm has a com- 
plexity scaling as N 1 . 

IV. CONCLUSION 

We have considered the formulation of generalized 
boundary conditions in a form useful for local elec- 
trostatic simulation algorithms. We have shown how 
to treat metallic boundaries for which one imposes a 
constant potential at a surface. An additional Monte 
Carlo move which transfers charge between conductors 
performs the Legendre transformation from a constant- 
charge to a constant-potential ensemble. 

The local algorithm (like Maxwell's equations) natu- 
rally generates a dipolar contribution to the solution to 



the Poisson equation in a periodic simulation cell. By 
choosing an anisotropic integration over the q = com- 
ponents of the electric field we reduce the spurious inter- 
action between different copies of a planar system. 

Molecular dynamics implementations of both prob- 
lems are potentially possible^. For the first problem 
of constant-potential simulation one could alternate be- 
tween a symplectic integrator (e.g. velocity Verlet) for 
the particles and a Monte Carlo step for charge trans- 
fer between conductors, or introduce a kinetic degree of 
freedom for the charges on each conductor. In the second 
problem of slab geometries one would need to follow the 
evolution of the electric field at each time step through- 
out the simulation cell, loosing the possibility of multi- 
scale updates. Implementation of these ideas has been 
performed in Ref. Ilia . 

Work financed in part by Volkswagenstiftung. 
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